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ABSTRACT 

We have developed a method of solving for multi-species chemical reaction flows 
in non-equilibrium and self-consistently with the hydrodynamic equations in an 
expanding FLRW universe. The method is based on a backward differencing scheme 
for the required stability when solving stiff sets of equations and is designed to be 
efficient for three-dimensional calculations without sacrificing accuracy. In all, 28 
kinetic reactions are solved including both collisional and radiative processes for the 
following nine separate species: H, H + , He, He + , He ++ , H~, , H%, and e~ . The 
method identifies those reactions (involving H~ and H%) ocurring on the shortest 
time scales, decoupling them from the rest of the network and imposing equilibrium 
concentrations to good accuracy over typical cosmological dynamical times. Several 
tests of our code are presented, including radiative shock waves, cosmological sheets, 
conservation constraints, and fully three-dimensional simulations of CDM cosmological 
evolutions in which we compare our method to results obtained when the packaged 
routine LSODAR is substituted for our algorithms. 

1. Introduction 



Hydrodynamical and microphysical processes of baryonic matter play an important role in 
structure formation at the smaller sub-cluster scales. Indeed, microphysics can play a more 
important role than gravity, especially when the cooling time of the gas is much shorter than the 
dynamical or Hubble times. Protogalactic clouds, which may form, for example, at high redshifts 
in CDM models can collapse through cooling instabilities to form an early generation of stars. The 
feedback from primordial stars can change the physical state of the pre-galactic medium and thus 
have considerable influence over the subsequent formation of stars and galaxies and the general 
state of the intergalactic medium (Couchman & Rees 1986; Tegmark et al. 1996). Microphysical 
processes are also very important at the center-most regions of cosmological sheets or pancakes. 
Originally studied by Zel'dovich (1970) in the context of neutrino-dominated cosmologies, sheets 
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are ubiquitous features in nonlinear structure formation simulations of CDM-like models with gas, 
and manifest on a spectrum of length scales and formation epochs. Cooling processes occur on 
very short time scales at the center-most densest parts of the pancake structures where stars and 
galaxies can form from the fragmentation of the gas. (Bond et al. 1984; Shapiro & Struck-Marcell 
1985; Yuan, Centrella & Norman 1991; Anninos & Norman 1996). 

It is well known that when nonequilibrium atomic reactions are properly taken into account, 
the cooling time can be shorter than the hydrogen recombination time and that gas which cools 
to 1 eV will likely cool faster than it can recombine. The effect of this nonequilibrium cooling is 
to leave behind a greater residual of free electrons and ions, as compared to the equilibrium case. 
The free electrons can be captured by neutral hydrogen to form H~ that subsequently produce 
hydrogen molecules. If large concentrations of molecules can form, the cooling is dominated by the 
vibrational/rotational modes of molecular hydrogen which acts to efficiently cool the gas to about 
1CT 2 eV, thereby reducing the Jeans mass of the gas. Hydrogen molecules can therefore play a 
crucial role in the formation of stars as they provide the means for cloud fragments to collapse 
and dissipate their energy. 

However, the typically high computational requirements and technical difficulties needed to 
solve the chemical rate equations relevant for H2 production in hydrodynamic flows have forced 
previous authors to impose simplifying assumptions such as the steady state shock condition 
which reduces the problem to zero dimension. In this case, only the time development of the 
hydrodynamic and thermodynamic variables are solved (Izotov & Kolesnick 1984; MacLow & 
Shull 1986; Shapiro h Kang 1987; Kang & Shapiro 1992). All of these studies have consistently 
found that the mass fraction of H2 can reach 10~ 3 behind sufficiently strong shocks, which is 
adequate to cool the gas to temperatures of order 10 -2 eV well within a Hubble time. More 
recently Haiman, Thoul & Loeb (1995) have investigated the formation of low mass objects in 
one-dimensional numerical calculations. They confirm the importance of H2 in the collapse of 
spherically symmetric isolated objects at high redshifts. 

Although much insight has been gained about the chemical aspects of molecular hydrogen 
formation and cooling, it remains to incorporate chemical reaction flows in realistic cosmological 
models. This paper discusses a method that we have developed for solving the kinetic rate 
equations with multi-species chemistry in nonequilibrium and self-consistently with the 
hydrodynamic and ./V-body equations in an expanding FLRW universe. The method is based 
on a backward differencing formula (BDF) for the required stability when solving stiff sets of 
equations, and is designed for both accuracy but especially speed so that it may be used in 
three dimensional codes with a minimal strain on computational resources. In all, we solve 
for 28 kinetic reactions including collisional and radiative processes for nine different species: 
H, H + , He, He + , He ++ , H~ , H2, H2, and e~, which we track individually with their unique 
mass transport equations. We have also implemented a comprehensive model for the radiative 
cooling of the gas that includes atomic line excitation, recombination, collisional ionization, 
free-free transitions, molecular line excitations, and Compton scattering of the cosmic background 
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radiation by electrons. 

The set of hydrodynamic, A^-body, kinetic and cosmological equations that we solve are 
summarized in §2.. Section 3. discusses the BDF method and its integration into the hydrodynamic 
solver. Several tests of our code are presented in §4., including ID radiative shock waves, 2D 
cosmological sheets, and fully 3D simulations of CDM cosmological evolutions in which we 
compare the BDF method to results obtained when the packaged routine LSODAR (Hindmarsh 
1983; Petzold 1983) is substituted in its place. We provide concluding remarks in §5.. 

Finally we note that a companion paper has been written which discusses the chemical model 
and cooling functions in more detail (Abel et al. 1996). In that paper, we motivate the model and 
argue its comprehensive nature in the choice of reactions, insofar as the formation of hydrogen 
molecules in cosmological environments is concerned. We also provide more up-to-date and 
accurate fits to the different rate coefficients. For completeness, we tabulate the reaction list and 
cooling processes in Appendices A and B of this paper, but leave the rate coefficients to Abel et 
al. (1996). 



2. The Equations 



The hydro dynamical equations for mass, momentum and energy conservation in an expanding 
FRW universe with comoving coordinates are 

^ + V-(p b v b ) + S-p b = 0, (1) 
ot a 

d(p b v b)i ) -. a 1 dp p b dcj) 

— ar - + v • [{pbVb > b] + 5 a pbVb - = -jWi ~ *Wi' (2) 

— » — * a 
— + V-(ev b )+pV-v b + 3-(e + p) = T-E, (3) 

where p b , p and e are the baryonic density, pressure and specific internal energy defined in 
the proper reference frame, v b is the comoving peculiar baryonic velocity, <ft is the comoving 
gravitational potential that includes baryonic plus dark matter contributions, a = 1/(1 + z) is the 
cosmological scale factor, and E and V are the microphysical cooling and heating rates. 

The equations for collisionless dark matter in comoving coordinates are 

dx d _ 



dv<i o a^ 1 

— = -2-v d - - 
at a a 



= -2-v d -—V(f>. (5) 



The baryonic and dark matter components are coupled through Poisson's equation for the 
gravitational potential 

V 2 = 47rGa 2 (p- p), (6) 
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where p = p b + Pd is the total density and p = 3i^Qr2o/(8vrGa 3 ) is the proper background density 
of the universe. 



The cosmological scale factor a(t) is given by Einstein's equation 



da 

dt~ H ° 



1 



fijif (- - 1) + «aK - 1) + 1 



a 



1/2 



(7) 



where VLm = ^b + ^d is the density parameter including both baryonic and dark matter 
contributions, S7a = A/(3Hq) is the density parameter attributed to the cosmological constant A, 
and Hq is the present Hubble constant. 

In addition to the usual hydrodynamic equations (1) - (3), we must also solve equivalent 
mass conservation equations for the densities pi of each of the nine separate atomic and molecular 
species that we track 

^ + V • ( Pt v h ) + 3% t = ±Y.H k A T )PjPl ± E hPh (8) 

3 I 3 

where the signs of each term on the right-hand-side depend on whether the process creates 
or destroys the species pi. The kji(T) are rate coefficients for the two body reactions and are 
functions of the gas temperature T. Explicit analytic fits for these coefficients over a broad range 
of temperatures, and a general discussion of the relevant chemical reactions, can be found in Abel 
et al. (1996). In all, we include 28 rate coefficients, one for each of the chemical reactions shown in 
Appendix A. The Ij in equation (8) are integrals due to photoionizations and photodissociations 

f°° l(v) 
Ij = totTjM-^-dv, (9) 

where T{y) = J-{v)/Att is the intensity of the radiation field, J-{y) is the flux, <Jj(v) are the 
cross-sections for the photoionization and photodissociation processes, and voj are the frequency 
thresholds for the respective processes. We note that the nine equations represented by (8) are not 
all independent. The baryonic matter is composed of hydrogen and helium with a fixed primordial 
hydrogen mass fraction of Hence we have the following three conservation equations 

hydrogen nuclei: p H + p H + + p H - + p H + + p H2 = p b f H , (10) 

helium nuclei: p He + Pne+ + PHe++ = Pb 0- ~ fn), (H) 

111 

charge conservation: p H + - p H - + -p H + + jPn e + + ^PHe++ = m H n e , (12) 

where n e is the number density of free electrons and tuh the proton mass. 

To complete the set of equations (1) - (3), we must also specify the equation of state 
appropriate for an ideal gas 
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where 7 = 5/3 is the ratio of specific heats for the baryonic matter, ks is boltzmann's constant, T 
is the gas temperature, and are the number densities for each of the different species. We also 
need to provide the necessary cooling E and heating V functions to the right-hand-side of equation 
(3) 

E = E Comp + J2J2^d T )PiP^ ( 14 ) 

3 I 

r = (15) 

3 

where Ec omp is the Compton cooling (or heating) due to interactions of free electrons with 
the cosmic microwave background radiation, and e.ji{T) are the cooling rates from two-body 
interactions between species j and I. The Jj(y) are integrals due to photoionizing and 
photodissociating heating 

J» = £ ^aj(u)l(u) {hV du. (16) 

We include a total of fourteen processes in the cooling function and three processes for heating. 
The physical mechanisms and mathematical expressions for each process are given in Appendix B. 



3. NUMERICAL METHODS 



It is well known that the differential equations describing non-equilibrium atomic and 
molecular rate reactions can exhibit variations on extreme time scales. Characteristic creation 
and destruction scales can differ by many orders of magnitude among the different species and 
reactions. As a result, explicit schemes for integration can be unstable unless unreasonably 
small time steps (smaller than the shortest dynamical times in the reaction flow) are taken, 
which makes any multi-dimensional computation prohibitively expensive. For this reason implicit 
methods are preferred for stiff sets of equations. These methods generally involve a Newton's 
iterative procedure to achieve convergence, and for large dimensional Jacobian matrices these 
implicit methods can also be very time consuming. A number of packaged routines exist which 
are based on identifying the disparity in time scales among the species and switching between 
stiff and nonstiff solvers. An example of such a package is the Livermore solver for ordinary 
differential equations with automatic method switching for stiff and nonstiff problems LSODAR 
(Hindmarsh 1983; Petzold 1983). However, an implementation of this solver in multi-dimensions 
is extremely costly in computer time and an alternative numerical scheme is desirable for fully 
three-dimensional calculations where computational speed is crucial. 



3.1. General Framework 
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We use an operator and directional splitting of the hydrodynamic equations (1) - (3) and (8) 
to update the fourteen state variables pb, Vb, e and pi. Six basic steps are utilized. 

First, the source step accelerates the fluid velocity due to pressure gradients and gravity, and 
modifies the velocity and energy equations to account for artificial viscosity 

Pb% = -v-Q, (18) 



'dt 

de 
dt 



-Q : W, (19) 



where Q is a second rank tensor representing the artificial viscous stresses (Stone &; Norman 1992). 
A staggered mesh scheme is utilized whereby the scalar variables pb, Pi, e, <ft and the artificial 
viscosity are zone centered, while the velocities are located at the zone interfaces. The pressure, 
potential, and viscosity gradients are thus naturally aligned with the momentum terms. 

In the second cooling/heating step, the energy changes are computed from "pdv" work and 
radiative cooling and heating from microphysical processes 

de — 

-^ = -pV-v-E cool + T. (20) 

We discuss solving this equation further in the following subsection. 

The third expansion step updates all state variables from the terms arising from the expansion 
of the universe 

dpb Q a / 01 x 

- -3-P6, (21) 



dt 

dpi _ Q d 

dt 



3-fH, (22) 



a 



dpbV -a _ 

~W = " 5 a^' (23) 

% = -Ae + P). (24) 
at a 

The homogeneous nature of the expansion allows a simple solution /o|, +A * = (a*/a' +A *) 3 p^, 
although a more generalized procedure is required for the energy equation in which an effective 
adiabatic index must be defined to eliminate the pressure term in the case of ionized gas (Anninos 
& Norman 1994). 



The fourth or transport step solves the advection terms 

/ PbdV = -J p b v-dS, (25) 
f PidV = -JpiV-dS, (26) 



dt 
d_ 
di 
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— J p b VidV = -J p b ViV • dS, (27) 
J edV = -Jev- dS. (28) 



d_ 

dt 



Several different monotonic schemes have been implemented, including donor cell, van Leer, and 
piecewise parabolic advection. All results presented here use the second order van Leer method 
which has the best accuracy-to-efficiency performance. 

A fifth step evolves the densities of the separate species according to the chemistry of the 
collisional and radiative kinetic equations 

| = ±EE^w±E% (29) 

3 1 3 

The methods for solving these equations is the focus of the next subsection. 

The total baryonic density p& and the density of each individual species pi are updated 
independently in the expansion, transport and chemistry steps. We are thus able to monitor the 
accuracy of our methods by evaluating the constraint equations for hydrogen, helium and charge 
conservation, equations (10) - (12). Over the course of a typical 3D calculation of approximately 
one-thousand time steps, we find maximum errors, which are mostly concentrated at the shock 
fronts, to be of order 10 to 30%. However, for increased stability and accuracy, we introduce a 
sixth step in our scheme to enforce the constraint equations at every time step. An important 
point that one must consider in taking this approach is that the errors can be larger than the 
concentration of those species that are depleted. It is therefore necessary to modify only the 
concentrations of the dominant species. For a hot gas in which H and He are mostly ionized, the 
hydrogen and helium constraints should be solved for the ionized components H + and He ++ . In 
cold (< 1 eV) gas the neutral hydrogen and neutral helium concentrations must be adjusted to 
satisfy the constraints. In this way, we are guaranteed to make only small (typically much less 
than one percent) fractional changes to any of the species at each timestep. 



3.2. Solving the Kinetic Equations 

More specific details of the numerical methods used in the hydrodynamic updates (17) - (28) 
can be found elsewhere (Stone & Norman 1992, Anninos et al. 1994). Here we emphasize the 
solution to the chemistry step (29). 

Notice that equations (29) can be written schematically as 

= d(T, nj ) - Di(T, n 3 ) m, (30) 

where n« = pi/^Aimn) and Ai is the atomic mass number of the ith element and mn is the proton 
mass. The Cj are the collective source terms responsible for the creation of the ith species. The 
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second terms involving Di represent the destruction mechanisms for the ith species and are thus 
proportional to ri{. Equation (30) suggests a backward difference formula (BDF) can be used in 
which all source terms are evaluated at the advanced time step. Discretization of (30) yields 

Lower order backward differentiation methods when applied to problems of the form y = f(y,t) 
are stiffly stable. This rather restrictive stability property is highly desirable when solving 
sets of stiff equations (Oran & Boris 1987). We have tried other less stable methods including 
higher order multi-step predictor-corrector schemes, various Runge-Kutta and Adams-Bashforth 
algorithms, and a Newton's procedure to solve the backwards differenced linearized equations. All 
of these alternative schemes have either proven to be unstable, less accurate or more expensive 
computationally compared to the simple BDF method. 

The solver can be optimized further by noting that the intermediaries H~ and in the 
molecular hydrogen production processes have large rate coefficients and low concentrations. They 
are thus very sensitive to small changes in the more abundant species. On the other hand, the low 
concentrations of H~ and implies that they do not significantly influence the more abundant 
electron, hydrogen and helium concentrations. This suggests that the nine species can be grouped 
into two categories: fast and slow reacting. The fast reacting group, comprised of H~ and H^, 
can be decoupled from the slower network and treated independently since the kinetic time scales 
for these species are much shorter than the characteristic times of the other seven species and 
the cosmological or gravitational times. H~ and can thus be considered in equilibrium at all 
times, independent of the hydrodynamic state variables. 

The expressions for the equilibrium abundances of H~ and H2 can be reduced by recognizing 
that reaction (19), according to Appendix A, can be neglected as a small order correction to 
H~ , due to the low concentrations of both species H~ and H^. Neglecting reaction (19), the 
equilibrium abundance of H~ can be written independent of H£ 

n _ = k 7 n H n e 

H k 8 n H + kun e + ki 5 n H + {he + k 17 )n H + +^23' 

where the variables k{ are the rate coefficients with subscripts referring to the reaction number in 
Appendix A. Then given n H - , the equilibrium abundance of can be written with no additional 
assumptions as 

n M+ = ; ; ; ; ; • I 33 ) 

2 k w n H + ki S n e + h 9 n H - + fc 25 + k 2 e 

The separation into fast and slow reacting systems helps to further increase the accuracy and 
stability of the BDF method when applied only to the slower network over the longer characteristic 
time scales required by the Hubble, hydrodynamic Courant, and gravitational free-fall times in 
cosmological simulations. 
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Due to the intrinsic nonlinearity of equation (30), not all source terms can be evaluated at 
the advanced time levels. Significant errors (of order 20% as measured by the final fractional 
abundance of hydrogen molecules) can be introduced if the source terms Q and are evaluated 
at the current time level at which the species data is known. Improvements to this crude 
approximation can be made by sequentially updating each species in order, rather than updating 
all species simultaneously from the data at the past time step. For example, the order in which 
we solve the rate equations was determined finally through experimentation to be H, H + , He, 
He + , He ++ and e~, followed by the algebraic equilibrium expressions for H~ and , then in 
the end updating H2 also using the BDF scheme (31). The updated concentrations of the ith 
(and previous) species are used as source terms in the equation for n^+i. Further improvements in 
accuracy and stability can be made by mimicking more closely a fully BDF scheme by subcycling 
the rate solve over a single hydrodynamic Courant time step. The subcycle time steps are 
determined so that the maximum fractional change in the electron concentration is limited to 10% 
per timestep, ie. 5r e = en e /h e with e = 0.1. 

We note that this same subcylcing procedure can be used to update the energy in the 
cooling/heating step in equation (20). The equation of state (13) couples the energy and pressure 
through the gas temperature, and although we have tried a Newton-Raphson iterative procedure, 
we found it to converge very slowly or sometimes not at all because the cooling/heating rates 
are strongly nonlinear and non-monotonic functions of temperature (Anninos & Norman 1994). 
For this reason we solve equation (20) with an explicit method that subcycles the cooling source 
terms. The timesteps for each subcyle are determined as ee/\E — T\, where e = 0.1 as in the rate 
equation solve. This algorithm has been tested to be both fast and accurate. 

The robustness of our methods has been verified by switching the order in which the rate 
equations are solved relative to the other updates. We have also experimented with the sensitivity 
to the temperature (ie. time centered, retarded and advanced temperatures) used in updating the 
rate equations. The results are stable and unchanged under all these sorts of permutations. 

4. CODE TESTS 

The numerical scheme described in Section 3. has been implemented in two separate codes: 
ZEUS-2D and HERCULES. ZEUS-2D is a two dimensional Eulerian hydrodynamics code originally 
developed by Stone & Norman (1992) and modified for cosmology by Anninos & Norman (1994). 
HERCULES is a three-dimensional hydrodynamics code derived from a 3D version of ZEUS, but 
modified for cosmology and generalized to include a hierarchical system of nested grids by Anninos 
et al. (1994). Extensive tests of the two codes can be found in the references provided above. In 
this paper, we present tests only for the new additions to the codes, namely the multi-species 
chemistry and the non-equilibrium cooling. Due to the intrinsic complexity of such systems, 
the diversity of tests is rather limited. We consider two crucial and relevant (for cosmology) 
tests: radiative shock waves and cosmological sheets. In addition we have also developed an 



independent method of solution that is based on the well-tested packaged solver LSODAR. Fully 
three-dimensional calculations of large scale structure formation are presented, comparing results 
from our method to that of LSODAR. 



4.1. 



Radiative Shock Waves 



Chevalier & Imamura (1982) and Imamura, Wolff &; Durisen (1984) have demonstrated 
through both analytic and numerical calculations that the fundamental mode of oscillations in 
one dimensional radiative shocks with cooling rates oc p 2 T a are unstable if a < 0.4. Because the 
cooling rate is determined from the density of electrons and ions, a comparison to their published 
work provides an excellent test of both our cooling algorithm and the reaction network. 

The initial data for this test is characteristic of pre-shock flows expected from the collapse of 
Zel'dovich pancakes 



corresponding to a uniform flow of gas along the — x direction. Reflection boundary conditions are 
imposed at x = and we use 100 zones to resolve a spatial extent of L = 2.43 x 10~ 4 Mpc. We 
assume only bremsstrahlung cooling of the form in Appendix B, which has an exponent of a ~ 0.5. 

A shock wave forms at the wall (x = 0) and propagates outward at a velocity v s ~ Ui n /3. 
As the heated gas cools, the shock begins to lose pressure support and slows down. Because 
the cooling rate is proportional to n e n H +T l l 2 (ie. oc /^T 1 / 2 ), the cooled gas experiences an 
accelerated energy loss as it gets denser. Eventually the higher density gas nearest the wall loses 
pressure support and the shock collapses and re-establishes a new pressure equilibrium closer to 
the wall. The shock front then begins again to propagate away from the wall, repeating the cycle 
of oscillations. This behavior is shown in figure 1 where the shock position x s is plotted as a 
function of time. Our numerical results indicate that the fundamental mode is indeed stable and 
damped, consistent with the analytic results of Chevalier and Imamura (1982) and the numerical 
simulations of Imamura et al. (1984). 

The shock jump conditions 



provide a more quantitative check of our numerical results. For the choice of initial data (34), the 
jump conditions give p post = 1.9 x 10 -24 g cm' 3 and p pos t = 1-8 x 10~ 10 g cm' 1 s' 2 , in excellent 
agreement with our numerical results p pos t = 1.9 x 10 -24 g cm' 3 &ndp post = 1.7 x 10~ 10 g cm' 1 s' 2 . 



p = 4.72 x 10~ 25 g cm' 3 , 
e = 1.0 x 10~ 30 g cm" 1 s~ 2 , 
v x = —Ui n = —1.7 x 10 7 cm s 



-i 



(34) 




(35) 
(36) 



' 1 

12 3 4 

Time(10 15 s) 

Fig. 1. — The location of the shock front x s as a function of time for the one dimensional radiative 
shock test. The complete nine species network of reactions is solved using the BDF method and 
a single source of cooling, Bremsstrahlung. The cooling time, maximum shock position, period of 
oscillation, and values of postshock state variables are all consistent with the analytic and numerical 
results of Chevalier & Imamura (1982) and Imamura et al. (1984). 

We can estimate the maximum distance the shock front will travel as 

%max = T coo iV s = -r- V s , (37) 

E 

where t coo i is the cooling time and v s is the shock speed. Substituting the bremsstrahlung cooling 
formula at the temperature predicted by the jump conditions and assuming a fully ionized gas with 
n e = n H + + 2n He ++ and hydrogen mass fraction $h = PhI Pb = 0.76, gives x max ~ 6.5 x 10 20 cm. 
This is again consistent with our numerical result 6.7 x 10 20 cm. 

Chevalier and Imamura also characterize their linearized analytic solutions by the frequency 
of oscillations ui = 0.31 in units of Ui n /x s where x s is the average shock front position. Defining 
the period P of perturbations as the time for the shock to first collapse back to the wall, we find 
P = 3.9 x 10 14 s and loj = (2n / P)(x s /ui n ) = 0.33, again in good agreement. 



4.2. Cosmological Sheets 
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We use the Zel'dovich (1970) solution to set up a linearized single mode perturbation for the 
collapse of gas in one dimension 

1 + z c sin(%) 



1 + z k ' 
H (l + z c )y/l + z 

l + Zc 



sin(kq) 



1 - 



cos(kq) 



k 
-i 



(38) 
(39) 
(40) 



1 + z 

where x and v x are the comoving positions and velocities, p the proper density, q the unperturbed 
coordinate, k = 2ir/\, A the comoving wavelength of perturbation, Hq the present Hubble 
constant, and z c the redshift corresponding to the collapse time. Parameters for the calculations 
presented here are the following 



1.0, 







B 



0.04, 



H = 70 km s _1 Mpc _1 , 



5, A = 10 Mpc. (41) 



At the higher temperatures (S> 1 eV) characteristic of shock heated gas in high velocity 
pancake structures, the kinetic time scales are extremely short compared to the Hubble time. 
Neglecting the photoionization processes, collisional ionization equilibrium is then a good 
approximation to the following pairs of reactions 



neutral hydrogen: H + e — > H + + 2e , 

H+ + e" H + 7, (42) 

neutral and singly ionized helium: He + e~ — ► ffe + + 2e _ , 

Fe+ + ->■ i3e + 7, (43) 

negatively charged hydrogen: H + e~ — ► + 7, 

H~ + e~ -> H + 2e", (44) 

ionized molecular hydrogen: + — > iJ^ + 7, 

# 2 + + e - 2H, (45) 

hydrogen molecules: + H —>■ H 2 + e~ , 

H 2 + e~ -> 2H + e~. (46) 



For a fully ionized gas, we have 

Ph+ 
Pb 



= Jh, 



(47) 
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PHe++ = l-fn, (48) 



Pb 

ne n H + + 2n He ++ _ 1+Jh 

n n H + + 4n He ++ 2 



(49) 



where fu = 0.76 is the mass fraction of hydrogen. The corresponding equilibrium fractional 
abundances from reactions (42) - (46) can then be written as 

(50) 
(51) 
(52) 
(53) 

f2fI (54) 



n H 
n 


k 2 


nue 


h n He + 


n 


k 3 n 


n H - 
n 


- k? k2 Ih 






n 


ki 8 fa 1 + f H ' 


n H2 
n 


k 8 k 7 / k 2 \ 2 
ki2 k 14 \kiJ 



1 + fH 

which are functions only of the gas temperature. 

In figure 2 we show the mass fractions Pi/pb throughout the pancake structure at redshift 
z = 4.8. At this time, the shock front is located at a distance of approximately 3 x 10 -3 Mpc 
from the central plane at x = 0, and is propagating outward at an average comoving velocity 
of about 110 km/s (355 km s~ 1 relative to the infalling gas). Two distinct cooling layers 
form, as evidenced by the thick solid line representing the gas temperature. The gas between 
3 x 10~ 6 Mpc < x < 1 x 10 -5 Mpc first cools mostly by atomic processes to a temperature 
of about 1 eV. The second colder layer at x < 10~ 6 Mpc results from cooling by hydrogen 
molecules which form from the residual electrons leftover from the nonequilibrium cooling through 
the first plateau at 1 eV. (We refer the reader to Anninos &: Norman (1996) for further details 
concerning the chemistry, dynamics and radiative cooling of cosmological sheets.) Although 
nonequilibrium effects are important as the gas cools through 1 eV, collisional ionization 
equilibrium is a good approximation to the species concentrations in the hot gas between 
2 x 10~ 5 Mpc < x < 2 x 10~ 3 Mpc. The different symbols plotted across this region in figure 
2 represent the equilibrium abundances given by equations (50) - (54). Differences between the 
numerical results and the equilibrium estimates are less than 0.5% for all species except H 2 which 
differ by roughly 5%. However, the larger discrepancy in the H 2 mass fraction is not due to 
numerical errors, but to the inadequacy of equation (54) to fully describe the kinetics. A more 
accurate equilibrium ratio is derived by considering all the reactions involving H 2 

n H2 k 8 n H -n H + k 10 n H +n H + k 19 n H +n H - 
n knn H + + ki 2 n e + ki 3 nn 

Equation (55) agrees with our numerical results to within about 0.1%. 
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Distance from Midplane (Mpc) 

Fig. 2. — Mass fractions of the different chemical species at redshift z = 4.8 across the cosmological 
sheet which first collapses at z = 5. Two distinctive cooling plateaus are evident in the temperature 
profile (thick solid line). The gas first cools to 1 eV mostly from atomic hydrogen line cooling (the 
region between 3 x < x < 1 x 10~ 5 Mpc) then cools further (the region x < 10~ 6 Mpc) due to 
hydrogen molecules formed from the residual electron fraction left behind by nonequilibrium effects. 
However, collisional equilibrium is a good approximation to the hot gas between the shock and the 
outer cooling plateau. The symbols plotted in this hot region represent the analytic equilibrium 
abundances and differ from the numerical results by less than 0.5%. 

We have also run this same problem using the LSODAR routines in place of the BDF method 
to solve for all nine species in full non-equilibrium. A comparison of the two results at redshift 
z = 4.8 is shown in figure (3). The symbols represent the BDF calculation and the various line 
types are the LSODAR results for the different species. The two methods agree to within about 
0.5% throughout all the different pancake layers, hot and cold. Notice also the excellent agreement 
in the mass fractions of H~ and i?2~, which is further justification of the hybrid model (in which 
the fast reacting species are singled out over cosmological dynamical times to be in equilibrium). 
Finally we point out that the fractional abundance of hydrogen molecules that form in the central 
cooled gas is consistent with the steady state shock calculations of Shapiro <fe Kang (1987). 



4.3. 3D CDM evolution 
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Fig. 3. — A comparison of results between the BDF method (symbols) and the fully non-equilibrium 
LSODAR algorithms (lines) for the cosmological sheet of the previous figure, also shown at redshift 
z = 4.8. The two methods agree to within 0.5% throughout both the hot and cold pancake layers. 

A final test is presented for an actual three dimensional cosmological calculation in which 
we compare our method of solving the rate equations to that of LSODAR. Both accuracy and 
computational efficiency are stressed in this comparison. 

The simulation is performed for a flat (Q,q = 1) model universe with baryon mass fraction 
f&6 = 0.035 and Hubble constant Hq = lOO/i km s~ 1 Mpc~ 1 with h = 0.65. The baryonic matter is 
composed of hydrogen and helium in cosmic abundance with a hydrogen mass fraction fjj = 0.76. 
The initial data is the Harrison-Zel'dovich power spectrum modulated with a transfer function 
appropriate for cold dark matter (CDM) adiabatic fluctuations and normalized to a bias factor of 
b = 1.0. We begin the simulations at redshift z = 50 and evolve to the present time at z = 0. The 
computational box size is set to L = 128 kpc (comoving) resolved by 32 3 cells. Our calculation 
thus has a spatial grid resolution of 4 kpc and baryonic mass resolution of 2.6 x 10 2 M Q which 
is just marginally adequate to resolve cooling flows and the formation of hydrogen molecules. 
The computational demands of LSODAR implemented in three dimensions prohibits comparative 
calculations of much higher resolution. 

We have run two separate simulations with identical initial conditions and model parameters. 
The only difference in the two runs is the rate equation solver: in the one case we use our BDF 
method, in the other LSODAR. The results of comparison are shown in figures 4 to 6. 
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Figure 4 shows the cell distributions (defined by counting up the number of cells at a 
particular binned range of values) for three key variables that would be particularly sensitive to 
errors in the rate solver: (a) gas temperature which is constructed from the concentrations of 
all the species combined, and (b) molecular hydrogen mass fraction pHil Pi and electron number 
density fraction n e /n. All results are shown at redshift z ~ 3. The prominent peaks in the H2 
and electron fraction distributions correspond closely to values initialized at the start of the runs. 
Because most of the box volume is comprised of cosmic voids that do not undergo shock heating, 
the molecular hydrogen and electron fraction do not change significantly in most of the cells. 
The relative rms deviations between the BDF and LSODAR results are 0.018, 0.039 and 0.13 for 
the gas temperature, molecular hydrogen and electron fractions, respectively. We note, however, 
that deviations are due in part to the final output redshifts not being exactly the same in the 
two runs; the relative difference is about 0.01. The larger differences found here (as compared to 
the cosmological sheet calculations) can also be attributed to the coarse grid resolution and the 
resulting large Courant time step in the hydro dynamical calculation. The agreements improve 
with resolution. 

Figure 5 shows contour graphs of the fractional volume of those cells with a particular 
combination of temperature and baryonic density, electron number density fraction and molecular 
hydrogen mass fraction. Two sets of graphs are presented: The three plots in the first column 
are results from the BDF method, the second column are the LSODAR results. We note that 
the sharp boundaries in the electron contours at the level log 10 (n e /n) ~ —3.3 correspond to the 
initialized fraction at z = 50: n e /n = 1.2 x 10^ 5 (/iSlf,) —1 (Peebles 1993). The volume weighted 
distribution is concentrated at the initial value for H2 since there is no mechanism to efficiently 
create nor destroy molecules at the low temperatures of the voids, and slightly lower than the 
initial value for the electron fraction since the expansion (cosmological and gravitational) of the 
voids continues to cool and recombine the gas. 

Contours of the spatial distribution of n e /n and ph2/ 1 Pb are shown in figure 6. The data is 
projected (and averaged) along the z-axis at redshift z = 3. The first row is the electron fraction, 
the second molecular hydrogen. The first column are the BDF results, the second LSODAR. 
Notice that hydrogen molecules form preferentially within the high dense filamentary structures 
but mostly in the knot-like intersections of the filaments. These are the highest density regions 
where the gas cools rapidly and the electrons are depleted by recombination. Peaks in the H2 
concentration therefore correspond to valleys in the electron distribution. In comparing results 
from the two solvers, BDF versus LSODAR, we see the distributions in both figures 5 and 6 are 
basically the same. 

In addition to the accuracy of the BDF method, another very important point that should 
be stressed is the amount of computational time required to solve the rate equations. The BDF 
run takes about 1.2 CPU hours with the full reaction network on the NCSA Convex C3880, in 
contrast to the equivalent LSODAR calculation which takes about 16.7 hours to complete on the 
same machine. The speedup of the BDF method over LSODAR is roughly a factor of 14. 



5. Summary 



We have developed and tested a new scheme to solve a system of stiff kinetic equations 
appropriate for chemical reaction flows in cosmological structure formation. Twenty-eight 
chemical reactions for collisional and radiative processes are included in our model, which tracks 
nine separate atomic and molecular species: H, H + , He, He + , He ++ , H~ , H^i -^2> an d eT . 
The reaction network is solved in a self-consistent manner with the hydrodynamic, N-body and 
cosmological expansion equations, and the accuracy of the solver has been verified by performing 
a series of test-bed calculations that includes radiative shock waves, cosmological sheets and 
monitering the conservation constraints. We have also implemented a publicly available and 
well-tested solver called LSODAR in place of our scheme and made direct comparisons of the 
different results in one, two and three dimensions. We find our methods are both fast and accurate, 
making fully three-dimensional calculations of non-equilibrium cosmological reaction flows feasible. 

We have incorporated the species solver into two separate cosmological hydrodynamic 
codes: a two-dimensional ratioed grid code to model Zel'dovich pancakes, and a more general 
three-dimensional nested grid cosmological code HERCULES. Applications to date include 
investigations of star and galaxy formation in cosmological sheets (Anninos & Norman 1996), 
primordial star formation in CDM models (Zhang et al. 1996), and simulations of the Lya forest 
(Zhang, Anninos & Norman 1995; Charlton et al. 1996). In the future, we plan to extend this 
work and develop a more sophisticated treatment of radiation to account for self-shielding and to 
more accurately model the microphysics of optically thick gas. 
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APPENDIX A: Chemical Reactions 

The following is a list of all chemical reactions that we include in our calculations. Further 
discussions and justification of the completeness of this set of reactions can be found in Abel et al. 
(1996). There we also include explicit formulae for the different rate coefficients. 

Collisional Processes: 
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APPENDIX B: Cooling and Heating Processes 

The cooling rates E and photoionization cross sections (Jj{v) included in our calculations. 
We use units of ergs cm _3 s _1 for the rat es, cm? for t he cross sections, and degrees Kelvin for 
temperature T. Also, T n = T/10™, a, L = (v/vq^) — 1 with vq^ being the threshold frequencies of 
the ith species, and k\, k 3 and k§ are the rate coefficients for the ionizing chemical reactions (1), 
(3) and (5) listed in Appendix A. 

Collisional excitation cooling (Black 1981; Cen 1992) : 
7.50 x 1(T 19 (1 + v^r 1 exp (-118348/T) n e n H 
9.10 x 10~ 27 (1 + VUy 1 T-°- 1687 exp(-13179/T) n 2 e n H e 

5.54 x 10~ 17 (1 + v^r 1 T-°- 397 exp (-473638/T) n e n He+ 

Collisional ionization cooling (Shapiro Sz Kang 1987; Cen 1992) : 
2.18 x lO" 11 kin e n H 
3.94 x 10" 11 k 3 n e n He 
8.72 x 10~ n k 5 n e n He + 

5.01 x 10~ 27 (1 + V%r X T-°- 1687 exp(-55338/T) n 2 e n He + 

Recombination cooling (Black 1981; Spitzer 1978) : 
8.70 x 10- 27 T 1 ' 2 T 3 -°- 2 (1 + T^ 7 )- 1 n e n H+ 

1.55 x 10~ 26 T - 3647 n e n He+ 

1.24 x 10- 13 T- 1 - 5 [1 + 0.3exp(-94000/T)] exp (-470000/T) n e n He + 
3.48 x 10- 26 T 1 / 2 r 3 -°- 2 (1 +T0- 7 )" 1 n e n He++ 

Molecular hydrogen cooling (Lepp & Shull 1983). 

Bremsstrahlung cooling (Black 1981; Spitzer & Hart 1979): 

1.43 x 10- 27 T 1 / 2 [l.l + 0.34exp(-(5.5-log 10 T) 2 /3)] n e {n H+ + n He+ + n He++ ) 

Compton cooling or heating (Peebles 1971) : 
5.65 x 10- 36 (1 + z) A [T - 2.73(1 + z)] n e 

Photoionization heating crossections (Osterbrock 1974) : 

a H (v) = 6.30 x 10" 18 {a 2 H + 1)~ 4 [exp (4 - 4(tan _1 a H )/a H )\ [1 - exp (-27r/ajr)] -1 
o He (y) = 7.42 x 10- 18 [l.m(a 2 He + l)" 2 - 05 - 0.66(a 2 He + l)" 3 - 05 ] 

a He+ {u) = 1.58 x 10~ 18 {a 2 He+ + l)" 4 [exp (4 - 4(tan- 1 a He+ )/a He+ )} [1 - exp (-27r/a //e+ )]- 1 
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Fig. 4. — The cell number distribution is shown for three state variables at redshift z = 3. (a) 
Temperature: The upward and downward pointing triangles are results from the BDF and LSODAR 
solver respectively, (b) Mass fraction of H2 and the number density fraction of electrons: The up 
and down triangles are BDF and LSODAR results for n e /n, while the left and right triangles are 
BDF and LSODAR results for pH2/ 1 Pb- The rms differences are 0.018, 0.039, and 0.13 for the 
temperature, H2 mass fraction, and electron fraction respectively. The errors are due in part to 
the output redshifts not being the same in the two cases, and to the large Courant timestep arising 
from the coarse grid resolution (32 3 cells). 
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Temperature (K) 

Fig. 5. — Contour plots for the cell distribution of the logarithm (log 10 ) of baryon overdensity, 
electron fraction n e /n, and molecular hydrogen fraction pwil Pb versus the log of the gas temperature 
at redshift z = 3. Ten contour levels are shown: (1 x 1(T 5 , 5 x 10~ 5 , 1 x 10~ 4 , 5 x 1(T 4 , 1 x 1(T 3 , 
5 x 10~ 3 , 1 x 10~ 2 , 5 x 10~ 2 , 1 x 1CT 1 , 5 x 1CT 1 ), which represent the fraction of cells with the 
specified values of the two physical quantities. The three graphs in the left column are results using 
the BDF method. The right column are the LSODAR results. The two methods agree very well. 
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Fig. 6. — Contour plots of the spatial distribution of the molecular hydrogen mass and electron 
number density fractions at redshift z = 3. The lower left panel is PH2/ Pb from the BDF method, 
and the lower right from LSODAR. The upper left panel is n e /n from BDF and the upper right 
from LSODAR. The contours cover the entire plane of the computational box and result from 
projecting (and averaging) the data along the z-axis. Six contour levels are shown for molecular 
hydrogen, (-4.78, -4.75, -4.72, -4.69, -4.66, -4.63), and eight for the electron fraction (-3.40, -3.39, 
-3.38, -3.37, -3.36, -3.35, -3.34, -3.33). Again, in comparing the BDF with the LSODAR results, 
the two methods can be seen to agree remarkably well. 



